On the radial distribution function of a hard-sphere fluid 
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Two related approaches, one fairly recent [A. Trokhymchuk et al. J. Chem. Phys. 123, 024501 
(2005)] and the other one introduced fifteen years ago [S. B. Yuste and A. Santos, Phys. Rev. A 
43, 5418 (1991)], for the derivation of analytical forms of the radial distribution function of a fluid 
of hard spheres are compared. While they share similar starting philosophy, the first one involves 
the determination of eleven parameters while the second is a simple extension of the solution of the 
Percus-Yevick equation. It is found that the second approach has a better global accuracy and the 
further asset of counting already with a successful generalization to mixtures of hard spheres and 
other related systems. 



The fluid of hard spheres (HS) is a key model in the 
use and application of statistical mechanics to the study 
of the thermodynamic and structural properties of real 
fluids. This is due not only to the relative simplicity 
of the intermolecular potential but also to the availabil- 
ity of rather accurate (albeit semi-empirical) expressions 
for the corresponding equation of state (EOS)i as well 
as the exact solution of the Ornstein-Zernike equation 
with the Percus-Yevick (PY) closure^ that provides an 
approximate analytical expression for the radial distribu- 
tion function (rdf). Yet, interest in this system is very 
much alive since exact analytical expressions both for 
its thermodynamic and structural properties have defied 
derivation. A recent example of this interest of particular 
relevance to our purpose here is the analytical expression 
for the rdf proposed by Trokhymchuk et al. (TNJH). 3 

Fifteen years ago, two of us4 developed an analytical 
method to determine reliable expressions for the rdf g (r) 
of an HS fluid. Such method provides a reasonable exten- 
sion of the solution to the PY equation by imposing weak 
physical requirements. Apart from being relatively sim- 
ple, requiring as the only input the EOS of the HS fluid, it 
has the further asset of including by construction the full 
thermodynamic consistence between the virial and com- 
pressibility routes to the EOS. The method uses a ratio- 
nal function ansatz for a function related to the Laplace 
transform of rg (r) , and so it is referred to as the Rational 
Function Approximation (RFA) method. It has subse- 
quently been used to study the metastable region past the 
liquid-solid transition^ successfully adapted and gener- 
alized to deal with other systems (such as HS mixtures^ 
sticky hard spheres^ and square- well fluids^), used to 
derive other structural properties (such as the contact 
values of the derivatives of the rdf, the direct correlation 
function, and the bridge function^) and also in connec- 
tion with the perturbation theory of liquids^ 

Interestingly enough, there are several common fea- 
tures in the RFA method and in the TNJH approach 3 
to derive an approximate analytical form for g (r) . For 
instance, the underlying philosophy is the same as ev- 
idenced through the following text that we quote from 



Ref. y: "A third possibility towards obtaining a nonem- 
pirical analytic representation of g(r) is to assume a cer- 
tain theoretically justified functional form for g(r) (e.g., 
a result of the solution of the OZ equation) with ad- 
justable parameters and impose then certain constraints 
(e.g., thermodynamic consistency) to determine these pa- 
rameters." Such a statement is totally applicable to the 
RFA method. Furthermore, both use the asymptotic be- 
havior of the rdf for r — > oo, they enforce thermody- 
namic consistency, and both also require as input only 
a prescribed compressibility factor Z = p/pfc^T and its 
associated isothermal susceptibility x = [d(pZ)/dp)~ 1 . 
Therefore it seems natural to compare the outcomes of 
both approaches in order to assess their merits and limi- 
tations. The purpose of this Note is to carry out such a 
comparison. 

In the TNJH approach? 3 , the rdf is written as 
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where r m (1 < r m < 2) denotes the first minimum of g(r) 
and the "depletion" (d) and "structural" (s) parts have 
the forms 

5d(r) = -e^-D + £ cos (/3(r - 1) + 7 )) e"^ 1 ' (2) 
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respectively. In Eqs. CEJ-© an d henceforth units are 
taken in which the diameter a has a value of 1. The 
depletion and structural parts have the same functional 
dependence on r as the exact solution of the PY equa- 
tion in the shell 1 < r < 2 and in the limit r — > oo, 
respectively. The approximation Q contains eleven pa- 
rameters to be determined as functions of the density p 
(or, equivalently, of the packing fraction r\ = irpa 3 /6). 
Trokhymchuk et al. impose five natural conditions: con- 
sistency with prescribed Z and x, continuity of g(r) and 
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its first derivative at r = r m 
at r — r m . More explicitly, 



and the minimum condition 
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In Ref. 3 the Carnahan-Starling-Kolafa (CSK)U EOS is 
used for Z and x- Furthermore, the parameters 7 and p 
are fixed at their PY expressions^ k and u> being given 
by parameterizations obtained by Roth et alH This still 
leaves seven parameters to be determined under the five 
constraints J3J| and By using the first equation of 
(J2J and the first and third equations of (J5J), A, B, C, and 
5 are expressed in terms of a, [3, and r m , as well as of 
g m = g(r m ). These remaining four unknowns are deter- 
mined by applying two conditions, namely, the second 
equation of (@J and the second equation of ©, and "by 
minimizing an appropriate functional," the form of which 
is not given in Ref. 3. In any case, due to the highly non- 
linear character of the problem, the numerical solutions 
for a, /3, r m , and g m are parameterized as functions of r\ 
in the range 0.2 < pa 3 < 0.9. 

On a different vein, in the RFA method the Laplace 
transform of rg (r) is taken to be given by 



G(t) = C{rg(r)} 



t 
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where 



$ (t) = (1 + Sit + S 2 i 2 + S 3 t 3 + S 4 i 4 )/(1 + Lit + 
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L 2 t 2 ), 



the six coefficients Si, £2, 1S3, 54, Li, and Li (which de- 
pend on the packing fraction) being evaluated in an al- 
gebraic form by imposing the following requirements*^ 

(i) x must be finite and hence the first two integral mo- 
ments of the total correlation function h (r) = g (r) — 1, 
i.e., J^°drr n h(r) with n = 1,2, must be well defined; 

(ii) the approximation must be thermodynamically con- 
sistent with a prescribed EOS, i.e., Eq. (0} must be ver- 
ified. Using the first requirement one finds that Li, Si, 
£2, and £3 are linear functions of Li and S4. The choice 
Li = S4 = yields the PY solution^ On the other hand, 
imposing the requirement (ii) leads^ to explicit expres- 
sions for Li and S4 in terms of x an d Z. In turn, the rdf 
is given by 



9 (r) = Hv ^ Vn (r ~ n) 9 (r 
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with O(r) the Heaviside step function and (p n (r) = 
— [$ (t)] - "}, which can be explicitly expressed by 
using the residue theorem £ To close the problem one has 
to give an EOS to fix Z and x, so all the procedure is a 
function of this choice. Here, we take the same choice as 
in Ref.0 namely the CSK EOSii 



The form of the TNJH approximation, Eqs. 
is particularly simple, while the RFA method, as hap- 
pens with the PY solution, gives an explicit expression 
for each separate shell n < r < n + 1 . On the other 
hand, the formal simplicity of the TNJH is at the ex- 
pense of introducing a number of unknowns exceeding 
the number of constraints, so that minimization of a cer- 
tain empirical functional is needed. As a consequence, 
the problem becomes highly nonlinear and the numerical 
solutions are fitted to given expressions; a change in the 
desired EOS would then require to start all over again 
the fitting procedure. In contrast, the RFA only requires 
six coefficients which are explicitly expressed in terms of 
an arbitrary EOS. From a less practical but more fun- 
damental point of view, it can be added that the TNJH 
g(r) introduces an artificial second-order discontinuity at 
the minimum r = r m , while it forces g(r) to be analytic 
at r = 2, 3, . . . 



0.01 



a 0.00 

bo 



If -0.01 h 

0.02 



s 0.00 

60 



p=0.4 



-RFA 
TNJH" 



p=0.6 



p=0.9 




FIG. 1: (Color online) Relative deviation Ag(r)/ gMo{r) = 
[s( r ) — Sto(f)] /<7md(0 of the values of the rdf of an HS fluid 
with respect to the molecular dynamics simulation data,— 4 
jmdH, as a function of r for pa 3 = 0.4 (top panel), pa 3 = 0.6 
(middle panel), and pa 3 — 0.9 (bottom panel). The dashed 
lines represent the results of the approach of Ref. y, while the 
solid lines refer to those of the RFA methodi- 



We are now in a position to carry out the aforemen- 
tioned comparison. Since both methods lead to reason- 
ably accurate values for the rdf, in Fig. Q] we present 
the results of the relative deviations of such values with 
respect to the very recent simulation data— for three 
characteristic densities. It is quite clear that only in the 
case pa 3 = 0.9 and for small distances, the RFA values 
give a somewhat poorer agreement with the simulation 
results. In any case, and recognizing the value of the 
TNJH approach, 3 it is fair to say that the RFA method 
provides a rather accurate and simple alternative to com- 
pute the structural properties of the HS fluid. Certainly, 
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there are fewer parameters in the RFA method and they 
have a clear physical link. Moreover, as the authors of 
Ref.0point out, the availability of analytical results in in- 
stances where simulations may fail (such as HS mixtures 
of a large size ratio) is crucial. In this context, the fact 
that we already have a simple extension of our approach 
to deal successfully with these and related systems£i& 
indicates that the RFA method is a valuable tool for the 
theoretical investigation of the structural properties of 
hard-core systems. 
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